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Within the recently introduced auxiliary master equation approach it is possible to address steady 
state properties of strongly correlated impurity models, small molecules or clusters efficiently and 
with high accuracy. It is particularly suited for dynamical mean held theory in the nonequilibrium 
as well as in the equilibrium case. The method is based on the solution of an auxiliary open quantum 
system, which can be made quickly equivalent to the original impurity problem. In its hrst imple¬ 
mentation a Krylov space method was employed. Here, we aim at extending the capabilities of the 
approach by adopting matrix product states for the solution of the corresponding auxiliary quantum 
master equation. This allows for a drastic increase in accuracy and permits us to access the Kondo 
regime for large values of the interaction. In particular, we investigate the nonequilibrium steady 
state of a single impurity Anderson model and focus on the spectral properties for temperatures 
T below the Kondo temperature Tk and for small bias voltages (j>. For the two cases considered, 
with T « Tk/4: and T « TkI 10 we hnd a clear splitting of the Kondo resonance into a two-peak 
structure for 0 close above Tk- In the equilibrium case {(f = 0) and for T « TkI'I, the obtained 
spectral function essentially coincides with the one from numerical renormalization group. 

PACS numbers: 71.15.-m,71.27-l-a,72.15.Qm,73.21.La,73.63.Kv 


I. INTRODUCTION 

The equilibrium properties of the single impurity 
Anderson model (SIAM) and the associated Kondo 
model,A.2i originally devised in the process of investi¬ 
gating metal hosts with dilute magnetic impurities 
are nowadays well understoodRenormalization group 
(RG) methods provided first perturbative analyses;^ and 
especially the development of Wilson’s numerical RG 
(NRG)^ allowed to properly capture the universal low- 
energy physics, governed by an exponentially small en¬ 
ergy scale, the Kondo temperature Tjf 4 The field of cor¬ 
related impurity models has gained renewed interest due 
to novel experimental realizations in quantum dots,— 
single molecule transistorsand from a theoretical 
point of view, due to its importance for dynamical mean 
field theory (DMFT).— The extension of DMFT to 
the nonequilibrium case can be carried out within the 
Keldysh formalism— Nonequilibrium DMFT and dif¬ 
ferent applicable impurity solvers have been thor oug hly 
discussed in other work, see for instance Refs. [29I - I3 a 

In the present study we want to focus on the physics 
of the impurity problem out of equilibrium itself, with 
an implementation of the auxiliary master equation ap¬ 
proach (AMEA) ^^i^^ based on matrix product states 
(MPS). Already in a first study, where Krylov space 
methods were employed^ AMEA has proven to feature 
a systematically improvable accuracy and to yield a well- 
defined Kondo peak in equilibrium together with a split¬ 
ting in the nonequilibrium case. However, the exponen¬ 
tial scaling of Krylov space methods with system size sets 
a “hard limit” to the achievable accuracy, and thus to the 
lowest temperatures accessible. The MPS extension pre¬ 


sented here turns out to be crucial in order to achieve 
highly accurate results in the Kondo regime down to low 
temperatures and up to large interactions. In the equi¬ 
librium case, the accuracy of our results becomes even 
comparable to NRG. 

Specifically, we investigate the nonequilibrium steady 
state dynamics of a SIAM, which is driven by the 
coupling to two leads at different chemical potentials, 
caused by an external bias voltage (f). Impurity mod¬ 
els in such a setup were considered already by many 
groups, numerically as well as analytically— AS To 
give a brief non-exhaustive overview, different tech¬ 
niques employed are the noncrossing approximation,— AS 
real-time diagrammatic methods^Si Keldysh perturba¬ 
tion theoryjSS Keldysh effective field theory^iSS dual 
fermionsperturbative RG,— flow equations^i^ 
functional RGj^i^ real-time RG,—AS time-dependent 
density matrix RG^AS NRG^AS Monte Carlo meth- 
odS)S2AS as well as cluster approaches— The properties 
of the correlated impurity have been established in cer¬ 
tain limits, for example for high temperatures T Tk 
or high biases (f Tk, where the Kondo effect is 
strongly suppressed by decoherence and the problem re¬ 
duces to a weak coupling one.—iSSiSSiliAS A splitting of 
the Kondo peak in the spectral function was found at 
sufficiently high bias voltages and low T, with a two- 
peak structure pinned to the chemical potentials of the 
leads.ii.A5>^^^.68-70 (f <s: Tk and 

T <C Tk, linear response as well as Fermi liquid theory is 
applicable , 49,58,74 , 75 Nevertheless, the intermediate and 
low-energy nonequilibrium regime, where both T and cj) 
are of the order of and especially below Tk, remains chal¬ 
lenging and the spectral properties could not yet be com- 
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pletely resolved. Work in this direction has for example 
been done in Refs. [Hiiii. However, the extension of 
NRG to the nonequilibrium case still leaves open ques- 
tionSf^^ and the Monte Carlo approaches, even though 
numerically exact, are either limited to relatively high 
temperatures and short times, or involve a demanding 
double analytical continuation<22iZ^ With the work pre¬ 
sented here, we want to contribute to these findings and 
present well-resolved spectral data for cases where both 
T Tr- and (j) < Tr-. 

II. MODEL AND METHOD 

The basic idea of AMEA is to map a general correlated 
impurity model in or out of equilibrium, here referred 
to as the physical impurity model (IMph), onto an ap¬ 
propriately chosen auxiliary one (IMaux)) which is small 
enough to be solvable precisely by numerical techniques. 
The self-energy of IMaux serves then as an approxima¬ 
tion to the one of IMph. Specifically, IMaux is modeled 
by an open quantum system described by a Lindblad 
equation, which consists of a finite number of bath sites 
and additional Markovian environments. In the mapping 
procedure, the bath parameters of IMaux are optimized 
in order to reproduce the dynamics of IMph as closely 
as possible. By increasing the number of bath sites Nb , 
more optimization parameters are available and a conver¬ 
gence (typically exponential) towards the exact solution 
of IMph is achieved. The mapping is formulated in terms 
of the hybridization function of IMaux, which is obtained 
through a single-particle calculation, and the many-body 
problem is solved thereafter. 

AMEA itself and a solution strategy for the correlated 
IMaux based on exact diagonalization (ED) was presented 
in detail in Refs. Here, we make use of MRS in 

order to solve for the correlation functions, which enables 
us to treat auxiliary systems with a larger number of bath 
sites. In the following we briefly summarize the governing 
equations in AMEA and point out modifications in the 
construction of IMaux favorable for an MPS treatment. 
After that, the MPS implementation is discussed. 

A. Keldysh Green’s functions 

In general, for nonequilibrium situations Green’s func¬ 
tions are conveniently defined on the Keldysh con¬ 
tour Since we are particularly interested in the long¬ 
time limit, where a steady state is reached, time transla¬ 
tional invariance applies and the Keldysh Green’s func¬ 
tions can be written in the frequency domain 

_ ^G«(a.) G^(a;)\ 

k(w) - Q j , (1) 

with , and we denote by an underscore ^ a 

2x2 object in Keldysh space. Only in an equilibrium sit¬ 
uation the Keldysh component is related to the retarded 


one via the fluctuation dissipation theorem 

G^{uj) = 2i{l-2pM^,fi,T))^ {G^(c*^)} , (2) 

where pfd{uj, p,T) represents the Fermi-Dirac distribu¬ 
tion. In contrast, in a general nonequilibrium situation 
a distribution function is not known a priori and the 
Keldysh and the retarded component have to be con¬ 
sidered as independent functions. 

It is convenient to introduce the steady state lesser and 
greater Green’s functions 

G<{t) = t (ct(t)c) , G>{t) = -t (c(t)ct) , (3) 

for generic fermionic creation/annihilation operators 
/c, which are related to G^ and G^ by 

G^{io) - G^{w) = G>{uj) - G<{uj) = -2iTrA{uj ), 

G^(w) =G>(w) + G<(a;), (4) 

and A(uj) is the spectral function. Throughout this work 
we consider solely steady state expectation values and 
denote them in compact notation by (...), cf. Eq. ([ 5 ]). 

B. Physical impurity model 

In this work, we consider for IMph a single impu¬ 
rity Anderson model in a nonequilibrium setup, given 
by an impurity Hamiltonian two noninteracting 

fermionic leads representing the electronic reservoir Hres, 
and an impurity-reservoir coupling i/coup: 

I^ph — Hjnip 4“ I^res 4“ I^coup • (5) 

The correlated impurity consists of a single level with 
energy Sd and on-site Hubbard interaction U 

=ed Y. 4^. + U (d\d^ - 1) (dY - 1) , 

<^£{ 4 . 4 .} ^ ^ 

_ . ( 6 ) 

where d^^/d^ are fermionic creation and annihilation op¬ 
erators on the impurity site. The reservoir Hamiltonian 
can be written in terms of the energy levels exk and po¬ 
tentials £a for the two leads A 

Hres = Yj ^ ^Xka^Xka > (’^) 

Xe{L,R} k<7 

and the impurity-reservoir coupling is given by 

H^coup = ^ tx (^^Xka^a + h.C.^ , (8) 

with representing creation and annihilation 

operators for lead electrons. 

Throughout this work we consider the particle-hole 
symmetric case with Ed = 0, ELk = £Rk- 
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An externally applied bias voltage (j) results in an anti- 
symmetrical shift of the chemical potentials 
In Sec. IIII Cl we further consider for the on-site energies 
the case ei^/n = ±-1, whereas in the rest of the work the 
voltage does not shift the lead energies. This is irrelevant 
for (j) much smaller than the bandwidth. 

The Green’s function of IMph is given by the Dyson 
equation 

Gph (^) = ffo- ^ph(^) - ■ (9) 

Here, denotes the noninteracting Keldysh Green’s 

function of the decoupled impurity, i.e. = {lu — ed)~^ 
and {g~^)^ can be neglected. The hybridization function 
Apjj is given by the sum of contributions from the two 
leads 

Aph(w) = ’ (19) 

A 

where denote lead Green’s functions at the contact 
point in the decoupled case. Except in the calculations 
presented in Sec. IIII Cl we consider throughout this work 
a flat band model with the retarded component of g^ (w) 
given by 

-3KM} = ^0(D-|a;|), (11) 

where we choose the hybridization strength T = jD 
as unit of energy and take D = 10 T. The real part is de¬ 
termined via the Kramers Kronig relation. For the fit in 
the mapping procedure (see Sec. Ill D|) it is of advantage to 
deal with smooth functions of w, so that we introduce in 
Eq. pip a smearing of the cut-offs in the Heaviside func¬ 
tion, determined by Fermi functions pfd(w, ±11,0.5 T) 
with an artificial temperature 0.5 T. Since this modifica¬ 
tion is well outside the scale of the impurity energies, it 
does not affect the low energy physics. 

The decoupled leads are in equilibrium, so that the 
Keldysh component g^i^) of each lead is given by Eq. ([2]) 
with the corresponding chemical potential g,x. The tem¬ 
perature T is taken to be the same in both of the 
leads. Notice that the Keldysh component is the only 
T-dependent quantity and results for different T shown 
below differ only in the smearing of the Fermi edge in 
g^{uj). In particular, we are interested in temperatures 
close to and below the Kondo temperature Tk- As for 
other methods, the low-temperature regime is most chal¬ 
lenging (cf. Sec. Ill Dl and App.[B]). For a Hubbard inter¬ 
action of 17 = 6 r, as considered throughout the work, 
one finds for the flat band model Tk ~ 0.2 

The remaining unknown quantity in Eq. (|9l) is the self¬ 
energy Spjj(a;), which cannot be determined exactly since 
IMph is interacting and of infinite size. This is evaluated 
by means of the mapping to IMaux- 


C. Auxiliary impurity model 

For IMaux we take an open quantum system of finite 
size, embedded in Markovian environments and described 
by a Lindblad equation for the system density operator 
P- 

= ( 12 ) 

The Lindblad super-operator L = Ch + Cd consists of an 
unitary part ChP = —i[Ha.ux, p] and the dissipator Cd &s 
described below 

Additionally to the original impurity site we consider 
Nb bath sites arranged in a linear geometry. For conve¬ 
nience we choose Nb even and the impurity site at the 
center, specified by the index /. The Hamiltonian for 
IMaux is given by 

T^aux = 

iju 

Here n/o- = with the fermionic cre¬ 

ation/annihilation operators and the {Nb -f 1) x {Nb + 1) 
matrix E couples only nearest neighbor (n.n.) terms, i.e. 
it is tridiagonal in the chosen geometry. To end up with 
a noninteracting bath we allow at most for Lindblad op¬ 
erators that are linear in The dissipator is then 

given hy^ 

CdP = 2^rg) U,pcl - i 

ijcr ^ / 

+ 2 X! (claPCja -\{p^ • (14) 

ija ^ ^ 

Both matrices of coupling constants and F^^^ are 
symmetric and positive definitei^ 

A key aspect in AMEA is that the bath parameters in 
the Lindblad equation are not determined within conven¬ 
tional Born-Markov approximations;^^— but are only 
used as fit parameters to optimally reproduce Apjj(w) 
by Aaux(‘^)i see Sec. | HD | 

Once the parameters of IMaux are determined, the 
many-body problem is solved (cf. Sec. Ill El) in order to 
obtain the interacting Green’s function 

^ux(w) = ffo ^(w) - - Saux(w) ■ (15) 

At this point it is convenient to set ^j,(w) = ^„x(aj) = 
F(a;), so that we obtain from Eq. (0 a very accurate 
result for the Green’s function of IMph. In this way, the 
U = 0 limit is recovered exactly. 

D. Mapping procedure 

In order to have a faithful representation of the dy¬ 
namics of IMph by IMaux, we need to fulfill Aaux(‘^) ~ 
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Apjj(a;) as closely as possible. For local quantities and 
correlation functions on the impurity, the influence of 
the bath is completely determined by the hybridization 
function only, independently of the specific bath geome¬ 
try. Therefore, the mapping becomes exact in the limit 
4aux(w) = Aph(w). To achieve ~ Aph(w), we 

minimize the mean squared error between them as a func¬ 
tion of the bath parameters in the Lindblad equation, i.e. 
the matrices E, and 

It is important to stress that a single particle calcu¬ 
lation is sufficient to determine A,^^,j(a;), for which the 
Green’s functions rea d^^i^^ 

G^{uj) =[uj-E + i (f(i) + F(2))) , 

G^{oj) = 2zG«(u;) (f^^) - F^^)) G^(cc). (16) 

Here, the inversion and multiplications are carried out for 
matrices in the site indices. The hybridization function 
is given in terms of the elements with impurity index / 

A«» = l/5o»-l/G^y;(u;), 

Afu.(a;) = . (17) 


D. With increasing number of bath sites Nb we observe 
that sharper features can be resolved. Therefore, a maxi¬ 
mal considered value of As converts to a lower bound for 
the ratio of temperature T to bandwidth 2D which can 
be reproduced by Ag^^,;(a;). More details on the mapping 
procedure are given in Add. [Bland Ref. 


E. Many-body solution 

1. Super-fermion representation 

As introduced in Refs. 1^1^ and made use of in 
Refs. Issll^ . the Lindblad equation m can be recast 
into a standard operator problem when considering an 
augmented fermion Fock space with twice as many sites. 
We use the notation of Ref. [9^ to which we refer to 
for further details, in combination with a particle-hole 
transformation in the “tilde” spacei^ The so-called left- 
vacuum reads 

|/) = E \M) 0 |{^) . (18) 

{'fT'icr } 


A single evaluation of the hybridization function is at 
most of 0{Nb^) and thus not time consuming. However, 
for a large number of bath parameters (> 20) the multi¬ 
dimensional optimization problem may become demand¬ 
ing and appropriate methods are needed. In particular, 
a parallel tempering approach has proven to be effective, 
which is discussed in some more detail in App. fAl 

Beyond the requirement Ag^^,;(a;) ~ Apj,(a;), complete 
freedom exists in choosing a suitable auxiliary system. 
For the many-body solution with MPS it is convenient 
to allow for nearest neighbor terms in the Lindblad cou¬ 
plings only, i.e. to restrict not only E but also the ma¬ 
trices fG) and F^^^ to a tridiagonal form. In this way 
one ends up with a geometry where the impurity couples 
to a bath with n.n. terms only. As discussed below, the 
bipartite entanglement entropy of IMaux can be reduced 
when imposing further that F^^j* has nonzero terms only 
for bath sites in one of the chains, e.g. for i,j > /, and 
on the other side, i.e. for i,j < /. For the latter 
restriction we found that it affects the quality of the fit 
only in a minor way but significantly improves the appli¬ 
cability of MPS. 

It is important to note that the relevant energy scale 
for the mapping procedure is not F but the bandwidth 
2D. For a certain IMaux, one can adjust to different F- 
values by simply rescaling all terms in E with index /, 
i.e. the hoppings to the impurity site, without chang¬ 
ing other properties of Ag^a,;(a;).— On the other hand, 
one can rescale the whole hybridization function by mul¬ 
tiplying the matrices E, F^) and F^^^ by the desired 
factor. Therefore, the complexity of the mapping pro¬ 
cedure is dominated by the smallest w-scale compared 
to the largest one. For the flat band model, this essen¬ 
tially means that one has to regard T and (f) in units of 


The summation runs over all possible many-body basis 

states \{nia}) of the original system and \{nia}) speci¬ 
fies those in the tilde system with inverted occupation 
numbers. N {{uia}) = total number of 

particles in state \{nicr}). 

The left-vacuum maps the density operator p{t) onto 
the state vector \p{t)) = p{t) |/). Thermodynamic expec¬ 
tation values are determined in this framework by expres¬ 
sions of the form {0{t)) = (/| O \p{t)). When evaluating 
(£p) 11) for the Lindblad equation (fl^ . one finds 

±\pit))=L\p{t)) , (19) 

where the super-operator C is replaced by an ordinary 
non-Hermitian operator L. In vector notation 

..., cj^,..., , (20) 


with and fermionic operators in the orig¬ 

inal and in the tilde system, respectively, the Lindblad 
operator L is given by 


zL = E' 


( 2 ) 


E + in 2F 
_2f(i) E-in 


Cg. — 2 Tr {E iA) 


+ U E ^/'^ + ’ (21) 


where n = F^^) - F^) and A = F^) + F^^). 
Clearly, L conserves the total particle number per spin 
J2i + ftirr)■ The Steady state |poo) = Ihn \p{t)) as 

t—¥oo 

well as |/) are situated in the half-filled, spin-symmetric 
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FIG. 1: (Color online) The upper part of the figure shows 
a schematic drawing of the auxiliary system in the super¬ 
fermion representation, for Nb bath sites and with the impu¬ 
rity located at the central site / = Nb/2. The upper chain 
corresponds to original sites and the lower chain to the addi¬ 
tionally introduced “tilde” sites, see Sec. Ill EH In the chosen 
fit restriction, see Sec. Ill Dl the coupling terms of the Lind- 
blad operator L, Eq. represent a ladder geometry with 
cross links. T^^^ causes a directional hopping from the upper 
to the lower chain and for it is vice versa. Moreover, 

is nonzero only for i,j < f and T^^ j* only for i,j > f. In the 
lower part, in grayscale, we schematically depict the tensor 
network for TEBD (Sec. IIIE2I) . where L is decomposed in 
nearest neighbor terms which are applied in an alter¬ 

nating manner. 


sector. Steady state expectations values and correlation 
functions are calculated 

{A{t)B) = {I\ Ae^*B Ipoo) , for t > 0 . (22) 

For tridiagonal matrices E, and see Sec. 

Ill Dl the coupling terms in Eq. (|2T]) represent a ladder 
system as depicted in Fig.[T] Sites on the original and 
the tilde system with the same index (i = j) ov i = 
j ± 1 are coupled with rates F-and F-^^- by a directional 
hopping. The restriction of F>J to the left side j < / 

and for F^^j* to the right side i,j > f leads to the situation 
that a circular current flows through the system. In this 
geometry one finds the tendency that sites on the left are 
filled in the original system and empty in the tilde system, 
whereas for the right side it is vice versa. This limits the 
possible hopping processes inside the chains and is in 
favor of a small bipartite entanglement entropyi^ 


2. Matrix product states 

A large amount of literature exists on MPS in gen¬ 
eral and for Lindblad-type problems in particular 
Here we briefly state the governing equations for the well- 
known MPS methods made use of in this work. 

We combine sites with the same index i in the original 
and in the tilde system to one “MPS-site”, with a local 
Hilbert space dimension d = 16 (see also Fig.[T|). For 


the resulting one-dimensional chain of sites it is straight¬ 
forward to write down a MPS representationiiSS 

Nb / d \ 

\p) = EK®*}) = n E A" ■ ( 23 ) 

{si} i=0 \si = l / 

Here, \p) is a generic many-body state with coefficients 
C{s^}. and A®* represents MPS matrices for site i with 
local quantum numbers The mapping Eq. (IMl) is 

exact for matrices which are exponentially large in Nb- 
However, even for much smaller matrix dimensions x 

dNB/2 

a very accurate representation of \p) is possible 
in many cases. For the auxiliary systems considered in 
this work, see also Sec. IHI A[ y Ri 1000 is sufficient when 
making use of Abelian symmetries of the Lindblad oper¬ 
ator Eq. m- Concerning the positivity of p, one should 
note that the form of Eq. (l23l) does not ensure it per 
construction !^^ However, we did not encounter un¬ 
physical results even for very small values of y. 

In order to calculate observables, a MPS representation 
of |/) is needed. One finds that Eq. (fTSl) can be recast into 
a state with y = 1, i.e. a product state, in which |/) is 
maximally entangled between original and tilde sites for 
the same index i. This is analogous to a purification of 
the identity operator/ 

When rewriting the Lindblad operator Eq. (|2T]) with 
tridiagonal matrices E, F^^^ and F*^^^ in form of a ma¬ 
trix product operator, one has couplings of n.n. sites 
only. This enables us to use very efficient time evolu¬ 
tion techniques as for example the time evolving block 
decimation (TEBD)ji2i Here, a Trotter decomposition is 
used to split the full time evolution exp(LAt) into small 
parts exp(Ti_i_|_i At) for neighboring sites, and terms with 
even and odd i are applied in an alternating manner, see 
also Fig.[TJ In this work we use splitting methods accu¬ 
rate to second order in We found that reducing 

the time step to At = 0.01 F“^ for the steady state and 
to At = 0.05 F“^ for the Green’s functions is usually suf¬ 
ficient. 

To obtain the desired steady state correlation functions 
of IMaux, for example G<, we proceed as follows: 

1. Calculate the steady state |poo) by time evolution 
with TEBD. Successively smaller time steps At are 
used in order to eliminate the Trotter error. Static 
observables and L \poo) = 0 may serve as conver¬ 
gence criterionii2i 

2. Apply Cj:^ to Ipoo) and time evolve the excited state 

to get G<{tn) = i{I\c^f„e^^’^Cf^\poG) at discrete 
points in the time domain. 

3. Employ linear prediction on the data G^{tn) 
and thereafter a Fourier transformation to obtain 
G^(u!) in the frequency domain 
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FIG. 2: (Color online) Temporal evolution of the bipartite 
entanglement entropy S in a typical IMauxwith Nb = 12, 
representing the case (fi = IT. The system is in the steady 
state for t < 0 and c/o- is applied to \pao) at t = 0. We show 
S{t) where it is largest, namely for the innermost bonds at 
the impurity, as well as for the next ones to the outside. 



w/r 


FIG. 3: (Color online) Spectral function in equilibrium: plot¬ 
ted for different number of bath sites Nb and compared with 
reference data from NRG.— Results are for [7 = OF and flat 
band leads, Eq. dill) , with D = 10 F and T = 0.05 F. 


III. RESULTS 

Before focusing on the nonequilibrium physics of the 
single impurity Anderson model, we briefly discuss the 
bipartite entanglement entropy of the auxiliary impurity 
model, and a benchmark for the equilibrium case. After 
that the spectral properties as a function of bias voltage 
are presented for two different temperatures, one clearly 
below and one above the Kondo temperature Tk- Fur¬ 
thermore, the bias dependence of observables such as the 
current and the double occupancy is discussed. In the 
last part of this paper a different density of states in the 
leads is considered, which allows to better resolve the 
physics at low temperatures and low bias voltages. 


A. Entanglement scaling 

Matrix product states are an efficient representation 
of many-body states with a low bipartite entanglement 
entropy S. The required matrix dimension y at a certain 
bond (i,i -|- 1) scales exponentially with the entropy at 
this bond, From Hermitian systems it is known, 

that ground states of gapped, one-dimensional systems 
obey an area law and are thus well-suited for MPS. Also 
an evolution in imaginary time converges well, but the 
real time evolution of excited states may become prob¬ 
lematic due to a build-up of entanglement— £ For the 
auxiliary impurity model investigated here, the behavior 
appears to be opposite. In general, the steady state |poo) 
of IMaux does not fulfill an area law and instead an in¬ 
crease of maxi with increasing system size Nb is 

observed— i Despite this, the time evolution of excited 
states is unproblematic, likely because of the damping 
involved, and the long time limit can easily be reached. 

We observe that the optimized parameters in IMaux 
strongly depend on the number of bath sites and on the 
external, physical parameters {(j), T, ...). Therefore, it is 
difficult to infer a reliable quantitative entanglement scal¬ 
ing with Nb- Qualitatively we hnd that max^ in¬ 

creases moderately with Nb and slower than linear. The 
magnitude of the entanglement is considerably reduced 
by the restricted setup for IMaux described in Sec. IIIDi 
which has the tendency towards a hlled and an empty 
bath chain in the steady state. In this setup, takes 

on the largest value at the central bonds which connect 
to the impurity site and falls off quickly with distance 
from the center. 

Independent of the actual scaling, the increase of bi¬ 
partite entanglement with Nb has the consequence that 
one is limited to certain system sizes. In this work we 
consider up to Nb = 16 with y = 1000, which is feasible 
in a reasonable amount of time. Most likely, one would 
need higher values of y in order to treat even larger sys¬ 
tems precisely. We checked the reliability of the results 
presented below by increasing the matrix dimension to 
y = 1500 in several cases, for different values of (f) and T. 
Furthermore, the time evolution for the Green’s functions 
was validated by reducing the time step to At = 0.01 F”^. 
Overall, we found in the worst cases relative differences 
in A{u:) up to 0(10“^). These errors are small enough 
for our purposes, so that we focus in the following on the 
accuracy of the mapping procedure, i.e. versus Nb- 

To analyze the temporal evolution of S, a typical time- 
dependent case is shown in Fig.O Here, t < 0 indicates 
the steady state regime and at t = 0 an annihilation op¬ 
erator is applied to |poo) in order to calculate the lesser 
Green’s function. To estimate the relevant time scale, 
Q{G^(f)} (not plotted) drops from 0.5 at t = 0 to 10“^ 
at t = 3F“^. As one can see, Si^i+i changes at first 
rapidly but saturates then and oscillates in time around 
a constant value. Thus, it is unproblematic to resolve 
G^{t) even on very large time scales. This was further¬ 
more checked for small Nb with an exact diagonalization 
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FIG. 4: (Color online) Bias-dependent spectral function (left) and retarded self-energy (right) for T = 0.05 F. Solid lines 
correspond to calculations with Nb = 16 and dash-dotted lines to Nb = 14, but in many cases they cannot be distinguished. 
Bias voltage (j) is given in units of F. Results are for 17 = OF and flat band leads, Eq. dill) , with D = lOF. 
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FIG. 5: (Color online) Bias-dependent spectral function (left) and retarded self-energy (right) for T = 0.5 F. Calculations are 
performed with Nb = 10 and other parameters are the same as in Fig.|4l 





FIG. 6: (Golor online) Double occupancy (left) and transport properties (right) as a function of bias voltage (j). Gurrent j is 
depicted with solid lines and the differential conductance djidfj) with dash-dotted lines. The latter is calculated with three-point 
Lagrange polynomials, based on the data for j as marked in the plot. Results are shown for T = 0.5 F with Nb ~ 10, and for 
T = 0.05 F with Nb = 14. Other parameters are as in Fig.|4] The linear response and equilibrium values for (nf-fUfi) are from 

NRG.-2£ 



























































solution as referenceJ^ When inspecting the short-time 
behavior of an asymmetry is evident. This results 

from the application of an operator to the original system 
alone, without changing the tilde system. 


B. Spectral and transport properties 

Before focusing on the nonequilibrium physics, we 
briefly present results for the equilibrium situation (j) = 0 
in Fig. [31 Here, a quasi-exact solution is provided by 
means of NRG.— For T = 0.05 F Ri Tk/4: the system 
is well inside the Kondo regime and the peak height of 
H(0)7rF R! 0.9 almost fulfills the T = 0 Friedel sum rule 
(H(0)7rF = 1 for T 0 ) 45 , 124,125 Results obtained with 
AMEA are shown for different system sizes, with partic¬ 
ular focus on the low energy physics. Noticeable differ¬ 
ences are apparent for Nb = 8, but, upon increasing the 
number of bath sites quick convergence is observed and 
excellent agreement with the NRG data is found. This 
shows that AMEA, especially with MPS, is a very ac¬ 
curate impurity solver also in the equilibrium case for 
T > 0. 

Regarding the accuracy of the calculations, we can 
state that Nb = 12 is essentially sufficient to provide 
reliable spectral data in equilibrium for T = 0.05 F. How¬ 
ever, for the nonequilibrium situations considered in the 
following one has to take into account that the accuracy 
of the mapping procedure is to some degree dependent on 
(j). This is analyzed in detail in App. |B]and here we solely 
want to note that the low bias ^ 1 F as well as the 

higher bias regime 4> > 2T converge more rapidly than 
the intermediate values, for T = 0.05 F. For the larger 
values of T used below the calculations are even easier, as 
one can achieve a very good mapping R Apjj(w) 

already for less than Nb = 12. 

After this benchmark, we now study the steady state 
nonequilibrium spectral properties for two different tem¬ 
peratures, one below and one above the Kondo tempera¬ 
ture. In Fig.jTjresults are presented for T = 0.05 F and in 
Fig. IS] for T = 0.5 F. In the first case, it is apparent that 
small bias voltages (j) <T cause a decrease and smearing 
of the Kondo peak, whereas larger voltages result in a 
splitting, see also Refs. Idll - lisll^roSl - lTOl It is known that 
with increasing current, resonant spin-flip scattering is 
prevented due to decoherence. Despite this, distinct ex¬ 
citations are clearly visible even at rather high bias volt¬ 
ages and located approximately at the positions of the 
chemical potentials Hl/r = This can be attributed 
to intra-lead processes, which however are strongly sup¬ 
pressed. In Fig.|4| we present furthermore the retarded 
self-energy, which enables us to better locate at which (j)- 
value splitting sets in. An upper bound can be estimated 
by the value (j> = 0.5 F, where —3{E^(u;)} exhibits two 
minima. In Sec. IHI Cl we resolve the physics at low bias 
in some more detail. 

In Fig. [5] the same system is considered for T = 0.5 F. 
As expected, the features are much broader and the 


Kondo peak for </> = 0 is strongly suppressed—^ De¬ 
spite of this, one can still note splitting and weak ex¬ 
citations a.t IJ.B/R = ^2 rather high voltages (^ > SF. 
In — S{S'^(w)}, only the result for </> = 2r exhibits two 
slight minima. One can thus infer, that the temperature 
dominates the decoherence processes on the impurity in 
this case and excitations at ^J.L/R are further suppressed 
and strongly smeared out. 

For both temperatures T = 0.05 F and T = 0.5 F, we 
present two observables of interest, the double occupancy 
and the current, in Fig.lH) The latter is obtained from 
the standard Meir-Wingreen expression— In the 
current it is obvious that the temperature strongly in¬ 
fluences the low bias regime, as is expected from linear 
response considerations. Especially the differential con¬ 
ductance enables to resolve the low-bias physics and we 
find a typical Kondo behavior . At higher voltages 
^ > 2r, however, one observes for T = 0.05F a slight 
increase of dj/dcj) due to charge fluctuations. At even 
higher voltages (</> > SF) both temperatures result in 
a similar linear current-voltage characteristic since the 
two spectral functions nearly merge into each other. The 
double occupancy (rif^nfi) exhibits an interesting be¬ 
havior for T = 0.05 F < Tk- In this case, (u/t^/ 4 .) and 
thus the charge fluctuation as well, exhibits a minimum 
(at (() R 2r). It originates from two competing mecha¬ 
nisms evolving with increasing (j): On the one hand, the 
enlarged transport window, approximately given by the 
interval (—f, f), increases and on the other 

hand, the suppression of resonant spin-flip scattering has 
the opposite effect. Apparently, the latter dominates 
initially at low bias. We observe a similar behavior in 
the temperature dependence of the double occupancy 
in the equilibrium case. We find a minimum 
in at T R 0.5 F. Therefore, the impurity is at 

this value in the local moment regime. When applying 
a bias voltage in the case of T = 0.5 F, the double occu¬ 
pancy increases monotonically with (j), as can be seen in 
Fig.|6| One can therefore conclude that the Kondo effect 
and its suppression with increasing cj) has a significant 
effect on the double occupancy. However, the particu¬ 
lar position of the minimum is not related to Tk but 
essentially only determi ned by the energy scale F, as dis¬ 
cussed in Refs. [T 2 I 1 M The reliability of our results is 
corroborated by the close agreement for ^ = 0 with the 
equilibrium values obtained by NRG (marked by circles 
in Fig.lHj). 


C. Low bias spectrum 

In order to better resolve the low energy spectral prop¬ 
erties of the Anderson impurity model, we now consider 
briefly the case of a Lorentzian density of states in the 
leads. In particular we replace Eq. (HU with 

9\i^) = iuJ-ex + ij)~^ , 


(24) 
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FIG. 7: (Color online) Bias-dependent spectral function (left) and retarded self-energy (right) for U = GT, T = 0.02F and a 
Lorentzian density of states in the leads, see Eq. (HU. Solid lines correspond to calculations with Nb = 16 whereas dash-dotted 
lines to ones with Nb = 14. The bias voltage <f) is in units of F. The inset on the left is for Nb = 16 and 0 = 0. 


where Eb/r = if. This can be produced by a bath 
consisting of one site with on-site energy e\ which fur¬ 
ther connects to a wide band given by 7 . We take 
r = —3{A^(0)} again as unit of energy and choose 
7 = S/ttF together with 2t '^= T. The Keldysh com¬ 
ponent is given by Eq. dU with = ±f, as before. 


A Lorentzian density of states is particularly suited 
for AMEA, since the retarded part alone can be 

fitted exactly with a single bath site. This simplification 
does not apply to the Keldysh component with its Fermi 
edges. Still, one can expect that the mapping procedure 
is more accurate than for a flat density of states and 
indeed we find that we are able to reproduce Apjj(w) by 
Aaux(‘^) more precisely for the same Nr- As a result, 
we can reach lower T with the same system sizes. For 
details on the achievable accuracy we refer to App. [BJ 


In particular we investigate the case T = 0.02 T and 
U = 6T and focus on bias voltages close to Tr. In ad¬ 
dition to the lower temperature especially the smaller ef¬ 
fective hybridization strength at the position of the Hub¬ 
bard bands leads to an increased separation of Kondo and 
Hubbard features in the spectral function, and thus, to 
an improved resolution. This can be seen in the peaked 
structure of the inset in Fig. [71 The smaller temperature 
allows us to analyze the behavior for lower bias voltages 
down to 0 = 0.1 r Also for this setup we find a similar 
dependence of the spectrum as a function of voltage as 
before, only at a decreased energy scale. The self-energy 
in Fig. [7] indicates that a splitting is first perceptible at a 
bias of 0 ~ 0.2 — 0.3 T. From our data we can thus con¬ 
clude that for bias voltages just above the Kondo tem¬ 
perature, a clear splitting of the Kondo resonance into a 
simple two-peak structure occurs. 


IV. CONCLUSIONS 

In this work we presented an improved formulation of 
AMEA, introduced in Refs. obtained by employ¬ 

ing matrix product states for the solution of the auxiliary 
master equation in the interacting case. This allowed us 
to treat larger auxiliary systems with more optimization 
parameters for the mapping pro cedure, as compared to 
the ED based solver in Ref. This is crucial, since 
the accuracy in AMEA increases exponentially with the 
number of optimization parameters. As a result, we ob¬ 
tained well-converged spectral data and static observ¬ 
ables, whose accuracy for the equilibrium case were com¬ 
parable to NRG down to low temperatures and for large 
interactions. More specifically, in the calculations pre¬ 
sented here, we were able to investigate the steady state 
properties of the single impurity Anderson model as a 
function of bias voltage 0 and at temperatures T well 
below the Kondo temperature Tr. In the spectral func¬ 
tion we obtained a prominent Kondo peak for 0 = 0 and 
T Ri TrIA, which compared very well to an equilibrium 
NRG calculation, and a broadening and subsequent split¬ 
ting of the peak when considering 0 > 0. Also for the 
case of a Lorentzian density of states in the leads, which 
enabled us to lower the temperature to T r; Tr/1Q, we 
found no evidence of a different behavior than a simple 
splitting of the Kondo peak. In order to locate the value 
of 0 at which the peak starts to split, it was advantageous 
to inspect the retarded self energy. From this we con¬ 
cluded that two excitations become visible for bias volt¬ 
ages just above the Kondo temperature, at 0 R 1 — 2 Tr. 

For the many-body solution with MRS it was of advan¬ 
tage to adjust the geometry of the auxiliary system and 
possible modifications were discussed. As in other stud¬ 
ies of Lindblad problems with MPS, we found an increase 
of the bipartite entanglement entropy S with system size 
Nr’^ However, the increase was moderate and slower 
than linear, which made it possible to treat auxiliary 
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FIG. 8: (Color online) Hybridization function as ob¬ 

tained from minimizing the cost function Eq. (HU with oJc = 
15 r and W{uj) = 1, for the flat band model Eqs. (HOI) . Hll) 
with (f) = 2V and T = 0.05E. Results on the left are for 
Nb = 12 and on the right for Nb = 14. 


open systems up to Nb ~ 16 reliably and within a rather 
short computation time (a couple of days). The value 
Nb = 16 is by no means a “hard limit” and much larger 
systems are expected to be feasible, especially when in¬ 
cluding additionally non-Abelian symmetries 

In general, the present MPS extension of AMEA con¬ 
stitutes a versatile and very accurate impurity solver for 
both equilibrium and nonequilibrium steady state situ¬ 
ations. Compared to the ED based solver presented in 
Ref. [ 3 ^ the computation time is longer but the achievable 
accuracy is much higher. Therefore, the MPS impurity 
solver is especially suited for situations in which a high 
spectral resolution is needed and a detailed investigation 
of the underlying physics is desired. 
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Appendix A: Multidimensional optimization 

In order to achieve Ag^^,j(a;) ~ Apjj(w), we optimize the 
bath parameters E, and For this a suitable 

parametrization is chosen, which yields a unique set of 
matrices E, P^^^ and P*-^^ for every parameter vector x. 
The mean squared error is quantified by a cost function 

UJc 

C{xf = f ^{^ph(w) - A“,,^(w;a;)}2lV(a;)dw, 

ae{R,K}_i^ 

with a certain cut-off Wc and weighting IV(w), which we 
take to be constant in the present paper. 

A variety of strategies exists to find the optimal pa¬ 
rameter set ajopt which minimizes a cost function as 
stated above. In previous work, Ref. we employed 
a gradient-based method with a large number of random 
starting points. Such a deterministic minimization works 
well for rather small problems, but becomes inefficient in 
the higher-dimensional case dim(a;) > 20. It is then of 
great advantage to employ methods which are able to 
overcome local minima. Appropriate Monte Carlo (MG) 
sampling based methods are for instance simulated an¬ 
nealing, multicanonical simulations or parallel tempering 
(PT)J^-237 Especially a feedback-optimized version of 
the latter has proven to be usefu l for our purposes. For 
details we refer to Refs. Il35lll36l and in the following we 
outline only briefly the implementation as used in this 
work. 

In PT, also called replica exchange, one regards C(a;) as 
an artificial energy, defines a set of artificial inverse tem¬ 
peratures and samples for each temperature from the 
Boltzmann distribution p^{x) = exp (—C(a;)/3m)- 

A replica x^ is assigned to each Pm and updated 
through a Markov chain with the Metropolis-Bastings al- 
gorithm.d^ These MC-sweeps generate a sequence of a;™, 
Z = 1, 2,..., which are distributed according to p^{x). In 
addition, a swapping of replicas a;)" and for neigh¬ 

boring inverse temperatures Pm and Pm+i is proposed 
after a certain number of sweeps. Again, a Metropolis 
probability is used for the swaps 

^ ^ (A2) 

with AC[" = (C(a;)") — C{xY^~^^)). The set of Pm in PT 
has the purpose that the low temperatures enable an effi¬ 
cient sampling of regions where C {x) is small and the ex¬ 
change with higher temperatures avoids trapping in local 
minima. To allow for an expedient exchange of replicas, 
the set of Pm needs to be adjusted. For our purposes we 
chose a feedback strategy which shifts the values Pm in 
order to achieve that the swapping probability Eq. (IA2I) 
becomes constant with respect to m. This strategy may 
not be the best possible choice in general, cf. Ref. Il37l 
but enables a fast feedback and quickly adjusts to large 
changes in the values C(a;)"). In addition, we modihed 
max(g™’'"~''^, gth) with a certain threshold 
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FIG. 9: (Color online) Convergence of with increasing Nb- Results with T = 0.5 F and T = 0.05 F are for the flat 

band case Eqs. (0,1111), and the ones with T = 0.02 F are for the Lorentzian density of states Eq. (1241) . For the cost function 
C, Eq.dHJ, we chose = 1 as well as cjc = 15r for the flat band model and cUc = 5r for the Lorentzian case. The 

normalization Co refers to the value of C for = 0. 



‘,'/r a,/r uiiT u}lT 


FIG. 10: (Color online) Convergence of the spectral function as depicted in Fig. [3] with increasing Nb, i.e. for the flat band 
model with [/ = 6 F and T = 0.05 F. 


probability {qth ~ 0.1), to avoid that during a PT-run 
a separation into several temperature sets occurs, which 
do not exchange replicas efficiently. This may violate 
balance conditions for thermodynamic observables but 
does not affect the applicability to minimization prob¬ 
lems. For the other PT-parameters we proceeded in the 
following way: In a single sweep each coordinate of a:™ 
was updated once and 10 sweeps were performed before 
attempting a swap. Around 20 — 30 inverse temperatures 
Pm were used. 


In general, one cannot expect to find the optimal solu¬ 
tion in a nontrivial high-dimensional problem, but with 
the PT algorithm as outlined above we obtain an Xmin 
which minimizes the cost function locally and may fur¬ 
thermore fulfill C(a:niin) ~ C(xopt) to good approxima¬ 
tion. For the largest systems considered in this work, 
Nb > 14, a good starting point was found to be impor¬ 
tant. For the case of tridiagonal E, and a 

convenient choice is to make use of Xmm from the next 
smaller system with Nb — 2. 


Appendix B: Convergence as a fnnction of Nb 

Figure [8] depicts two typical results of the optimization 
described in Ann.lAl for Nb = 12 and Nb = 14. It is ap¬ 
parent that rapid convergence is achieved when increas¬ 
ing Nb- For low temperatures T we find that the biggest 
error in occurs in the retarded component at 

the positions of the chemical potentials ^jll/r = 4=^, see 
Nb = 12. This is a consequence of optimizing 
and simultaneously. For higher temperatures, 

for instance T = 0.5 F, this effect is much less pro¬ 
nounced. 

A brief analysis of the convergence behavior of the 
mapping procedure with increasing Nb is given in Fig. 1^1 
We present values of the cost function C, Eq. (lAll) . for 
different temperatures and bias voltages. In general one 
finds an exponential convergence C oc exp(—rAs) to 
good approximation and the higher the temperature, the 
higher the rate of convergence r. By averaging over re¬ 
sults for different p we estimate a scaling of r oc Ti. One 
can deduce from the order of magnitude of C, that the 
calculations presented for A{u}) at T = 0.02 F (Fig. [ 7 ]) 
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are not converged to the same accuracy as the ones at 
T = 0.05 r (Fig.|31) or r = 0.5r (Fig.[S|), and larger sys¬ 
tems with Nb ^ 20 would be needed. However, the accu¬ 
racy is comparable to the Nb = 12 results for T = 0.05 F, 
which already yielded qualitative correct behavior and 
quite accurate spectral data, see also Fig. [3] and Fig.lTOl 
The influence of (j) is nonmonotonic and strongly depen¬ 
dent on the particular density of states in the leads. For 
the situations considered in this work we find the ten¬ 
dency that larger ^ result in larger values of C. For a 
more detailed analysis of the scaling with temperature 
and the mapping procedure in general we refer to Ref. 1^ . 

For the flat band case with T = 0.05 F we present 
a more thorough investigation by comparing the spec¬ 
tral function in the interacting case 17 = 6 F for different 
numbers of bath sites in Fig. [TUI As can be anticipated 
from the cost function C in Fig.jUl the cases cj) = 1/3 F 
and (j) = 2T are well-converged for Nb = 16, which man¬ 
ifests itself also in A{uj). The case cj) = T exhibits larger 
values of C and one can note more significant changes in 


Ajuj). Interestingly, rather high values of C are obtained 
for (/ = 4r, but nevertheless, the spectral function con¬ 
verges nicely. As discussed above for Fig.jHl the largest 
errors in correspond to short-scaled oscillations 

in These errors are likely to be averaged out 

once the spectral function exhibits rather broad features. 
This is exactly the case for higher bias voltages where the 
Kondo effect is strongly suppressed. On the whole, when 
inspecting Fig.lTOland also Fig.jSl one can note a slightly 
nonsmooth convergence with Nb, especially close to the 
Kondo regime for low (j>. This can be accounted to abrupt 
changes of spectral weight in around w = 0, 

when changing Nb- One possibility to suppress this ef¬ 
fect is to adjust the weighting function W(uj) in Eq. (jAll) 
accordingly. However, this is most probably only of im¬ 
portance when aiming to achieve even higher accuracies 
in A{uj) and with the choice Wjoj) = 1, one can regard 
the calculations presented in this work as unbiased and 
accurate over the whole w-domain. 
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